embryo_L1_L3_intestine_RNAseq
Rtpw
2/8/2022
Install libraries
Load libraries
source(file = "../../elt2_paper_functions.R")Analysis outline
- Import counts and format for DESeq2 analysis
- Filter low abundance counts based on counts per million (CPM)
- Filter for protein coding genes
- Perform differential expression analysis
- Visualize pairwise sample correlation
- Categorize gene expression with alternative hypothesis
- Visualize unshrunken pairwise differential expression
- Visualize shrunken pairwise differential expression
- Measure tissue expressed genes
- Visualize developmentally dynamic intestine genes
- Identify novel intestine genes in each stage for smFISH experiments
Import and format count data
# import counts
countsData <- read.delim(file = "../01_input/all.counts", sep = " ")
# preview counts
head(countsData)## chr start stop strand length embryo_cells_rep1
## WBGene00014450 MtDNA 1 55 + 55 0
## WBGene00014451 MtDNA 58 111 + 54 0
## WBGene00010957 MtDNA 113 549 + 437 0
## WBGene00010958 MtDNA 549 783 + 235 0
## WBGene00014452 MtDNA 785 840 + 56 0
## WBGene00014453 MtDNA 842 896 + 55 0
## embryo_cells_rep2 embryo_GFPminus_rep1 embryo_GFPminus_rep2
## WBGene00014450 0 0 0
## WBGene00014451 0 0 0
## WBGene00010957 0 0 0
## WBGene00010958 0 0 0
## WBGene00014452 0 0 0
## WBGene00014453 0 0 0
## embryo_GFPminus_rep3 embryo_GFPplus_rep1 embryo_GFPplus_rep2
## WBGene00014450 0 0 0
## WBGene00014451 0 0 0
## WBGene00010957 0 0 0
## WBGene00010958 0 0 0
## WBGene00014452 0 0 0
## WBGene00014453 0 0 0
## embryo_GFPplus_rep3 embryo_whole_rep2 embryo_whole_rep3
## WBGene00014450 0 0 0
## WBGene00014451 0 0 0
## WBGene00010957 0 0 0
## WBGene00010958 0 0 0
## WBGene00014452 0 0 0
## WBGene00014453 0 0 0
## L1_cells_rep1 L1_cells_rep2 L1_cells_rep3 L1_GFPminus_rep1
## WBGene00014450 0 0 0 0
## WBGene00014451 0 0 0 0
## WBGene00010957 0 0 0 0
## WBGene00010958 0 0 0 0
## WBGene00014452 0 0 0 0
## WBGene00014453 0 0 0 0
## L1_GFPminus_rep2 L1_GFPminus_rep3 L1_GFPplus_rep1
## WBGene00014450 0 0 0
## WBGene00014451 0 0 0
## WBGene00010957 0 0 0
## WBGene00010958 0 0 0
## WBGene00014452 0 0 0
## WBGene00014453 0 0 0
## L1_GFPplus_rep2 L1_GFPplus_rep3 L1_whole_rep1 L1_whole_rep2
## WBGene00014450 0 0 0 0
## WBGene00014451 0 0 0 0
## WBGene00010957 0 0 0 0
## WBGene00010958 0 0 0 0
## WBGene00014452 0 0 0 0
## WBGene00014453 0 0 0 0
## L1_whole_rep3 L3_cells_rep1 L3_cells_rep2 L3_cells_rep3
## WBGene00014450 0 0 0 0
## WBGene00014451 0 0 0 0
## WBGene00010957 0 0 0 0
## WBGene00010958 0 0 0 0
## WBGene00014452 0 0 0 0
## WBGene00014453 0 0 0 0
## L3_GFPminus_rep1 L3_GFPplus_rep2 L3_GFPminus_rep3
## WBGene00014450 0 0 0
## WBGene00014451 0 0 0
## WBGene00010957 0 0 0
## WBGene00010958 0 0 0
## WBGene00014452 0 0 0
## WBGene00014453 0 0 0
## L3_GFPplus_rep1 L3_GFPminus_rep2 L3_GFPplus_rep3 L3_whole_rep1
## WBGene00014450 0 0 0 0
## WBGene00014451 0 0 0 0
## WBGene00010957 0 0 0 0
## WBGene00010958 0 0 0 0
## WBGene00014452 0 0 0 0
## WBGene00014453 0 0 0 0
## L3_whole_rep2 L3_whole_rep3
## WBGene00014450 0 0
## WBGene00014451 0 0
## WBGene00010957 0 0
## WBGene00010958 0 0
## WBGene00014452 0 0
## WBGene00014453 0 0# print samples
colnames(countsData[6:ncol(countsData)])## [1] "embryo_cells_rep1" "embryo_cells_rep2" "embryo_GFPminus_rep1"
## [4] "embryo_GFPminus_rep2" "embryo_GFPminus_rep3" "embryo_GFPplus_rep1"
## [7] "embryo_GFPplus_rep2" "embryo_GFPplus_rep3" "embryo_whole_rep2"
## [10] "embryo_whole_rep3" "L1_cells_rep1" "L1_cells_rep2"
## [13] "L1_cells_rep3" "L1_GFPminus_rep1" "L1_GFPminus_rep2"
## [16] "L1_GFPminus_rep3" "L1_GFPplus_rep1" "L1_GFPplus_rep2"
## [19] "L1_GFPplus_rep3" "L1_whole_rep1" "L1_whole_rep2"
## [22] "L1_whole_rep3" "L3_cells_rep1" "L3_cells_rep2"
## [25] "L3_cells_rep3" "L3_GFPminus_rep1" "L3_GFPplus_rep2"
## [28] "L3_GFPminus_rep3" "L3_GFPplus_rep1" "L3_GFPminus_rep2"
## [31] "L3_GFPplus_rep3" "L3_whole_rep1" "L3_whole_rep2"
## [34] "L3_whole_rep3"# import metadata and process file
metadata1 <- read.table(file = "../01_input/RWP27_metadata.tsv", header = FALSE, stringsAsFactors = FALSE) %>% bind_rows(read.table(file = "../01_input/RWP26_metadata.tsv", header = FALSE, stringsAsFactors = FALSE)) %>%
bind_rows(read.table(file = "../01_input/RWP30_metadata.tsv", header = FALSE, stringsAsFactors = FALSE))
colnames(metadata1) <- c("Filename.Fwd", "Filename.Rev", "names")
head(metadata1)## Filename.Fwd Filename.Rev names
## 1 RW57_S10_L003_R1_001 RW57_S10_L003_R2_001 embryo_cells_rep1
## 2 RW58_S11_L003_R1_001 RW58_S11_L003_R2_001 embryo_GFPplus_rep1
## 3 RW59_S12_L003_R1_001 RW59_S12_L003_R2_001 embryo_GFPminus_rep1
## 4 RW60_S13_L003_R1_001 RW60_S13_L003_R2_001 embryo_whole_rep2
## 5 RW61_S14_L003_R1_001 RW61_S14_L003_R2_001 embryo_cells_rep2
## 6 RW62_S15_L003_R1_001 RW62_S15_L003_R2_001 embryo_GFPplus_rep2# separate and process sample info
metadata1 <- metadata1 %>% separate(names, sep = "_", into = c("stage", "sample", "rep"), remove = FALSE)
metadata1 <- metadata1 %>% mutate(stage = fct_relevel(stage, c("embryo", "L1", "L3")),
sample = fct_relevel(sample, c("whole", "cells", "GFPplus", "GFPminus")),
rep = fct_relevel(rep, c("rep1", "rep2", "rep3")),
names = fct_relevel(names, metadata1$names)
)
# Order columns according to metadata1 order
countsData <- countsData %>% select(chr:length, sort(metadata1$names))
head(countsData)## chr start stop strand length embryo_cells_rep1
## WBGene00014450 MtDNA 1 55 + 55 0
## WBGene00014451 MtDNA 58 111 + 54 0
## WBGene00010957 MtDNA 113 549 + 437 0
## WBGene00010958 MtDNA 549 783 + 235 0
## WBGene00014452 MtDNA 785 840 + 56 0
## WBGene00014453 MtDNA 842 896 + 55 0
## embryo_GFPplus_rep1 embryo_GFPminus_rep1 embryo_whole_rep2
## WBGene00014450 0 0 0
## WBGene00014451 0 0 0
## WBGene00010957 0 0 0
## WBGene00010958 0 0 0
## WBGene00014452 0 0 0
## WBGene00014453 0 0 0
## embryo_cells_rep2 embryo_GFPplus_rep2 embryo_GFPminus_rep2
## WBGene00014450 0 0 0
## WBGene00014451 0 0 0
## WBGene00010957 0 0 0
## WBGene00010958 0 0 0
## WBGene00014452 0 0 0
## WBGene00014453 0 0 0
## embryo_whole_rep3 embryo_GFPplus_rep3 embryo_GFPminus_rep3
## WBGene00014450 0 0 0
## WBGene00014451 0 0 0
## WBGene00010957 0 0 0
## WBGene00010958 0 0 0
## WBGene00014452 0 0 0
## WBGene00014453 0 0 0
## L1_whole_rep1 L1_cells_rep1 L1_GFPplus_rep1 L1_GFPminus_rep1
## WBGene00014450 0 0 0 0
## WBGene00014451 0 0 0 0
## WBGene00010957 0 0 0 0
## WBGene00010958 0 0 0 0
## WBGene00014452 0 0 0 0
## WBGene00014453 0 0 0 0
## L1_whole_rep2 L1_cells_rep2 L1_GFPplus_rep2 L1_GFPminus_rep2
## WBGene00014450 0 0 0 0
## WBGene00014451 0 0 0 0
## WBGene00010957 0 0 0 0
## WBGene00010958 0 0 0 0
## WBGene00014452 0 0 0 0
## WBGene00014453 0 0 0 0
## L1_whole_rep3 L1_cells_rep3 L1_GFPplus_rep3 L1_GFPminus_rep3
## WBGene00014450 0 0 0 0
## WBGene00014451 0 0 0 0
## WBGene00010957 0 0 0 0
## WBGene00010958 0 0 0 0
## WBGene00014452 0 0 0 0
## WBGene00014453 0 0 0 0
## L3_whole_rep1 L3_cells_rep1 L3_GFPplus_rep1 L3_GFPminus_rep1
## WBGene00014450 0 0 0 0
## WBGene00014451 0 0 0 0
## WBGene00010957 0 0 0 0
## WBGene00010958 0 0 0 0
## WBGene00014452 0 0 0 0
## WBGene00014453 0 0 0 0
## L3_whole_rep2 L3_cells_rep2 L3_GFPminus_rep2 L3_GFPplus_rep2
## WBGene00014450 0 0 0 0
## WBGene00014451 0 0 0 0
## WBGene00010957 0 0 0 0
## WBGene00010958 0 0 0 0
## WBGene00014452 0 0 0 0
## WBGene00014453 0 0 0 0
## L3_whole_rep3 L3_cells_rep3 L3_GFPplus_rep3 L3_GFPminus_rep3
## WBGene00014450 0 0 0 0
## WBGene00014451 0 0 0 0
## WBGene00010957 0 0 0 0
## WBGene00010958 0 0 0 0
## WBGene00014452 0 0 0 0
## WBGene00014453 0 0 0 0# Generate a table called "cts" out of the countsData table. Subset the countsData.
cts <- as.matrix(countsData %>% select(metadata1$names))
head(cts)## embryo_cells_rep1 embryo_GFPplus_rep1 embryo_GFPminus_rep1
## WBGene00014450 0 0 0
## WBGene00014451 0 0 0
## WBGene00010957 0 0 0
## WBGene00010958 0 0 0
## WBGene00014452 0 0 0
## WBGene00014453 0 0 0
## embryo_whole_rep2 embryo_cells_rep2 embryo_GFPplus_rep2
## WBGene00014450 0 0 0
## WBGene00014451 0 0 0
## WBGene00010957 0 0 0
## WBGene00010958 0 0 0
## WBGene00014452 0 0 0
## WBGene00014453 0 0 0
## embryo_GFPminus_rep2 embryo_whole_rep3 embryo_GFPplus_rep3
## WBGene00014450 0 0 0
## WBGene00014451 0 0 0
## WBGene00010957 0 0 0
## WBGene00010958 0 0 0
## WBGene00014452 0 0 0
## WBGene00014453 0 0 0
## embryo_GFPminus_rep3 L1_whole_rep1 L1_cells_rep1 L1_GFPplus_rep1
## WBGene00014450 0 0 0 0
## WBGene00014451 0 0 0 0
## WBGene00010957 0 0 0 0
## WBGene00010958 0 0 0 0
## WBGene00014452 0 0 0 0
## WBGene00014453 0 0 0 0
## L1_GFPminus_rep1 L1_whole_rep2 L1_cells_rep2 L1_GFPplus_rep2
## WBGene00014450 0 0 0 0
## WBGene00014451 0 0 0 0
## WBGene00010957 0 0 0 0
## WBGene00010958 0 0 0 0
## WBGene00014452 0 0 0 0
## WBGene00014453 0 0 0 0
## L1_GFPminus_rep2 L1_whole_rep3 L1_cells_rep3 L1_GFPplus_rep3
## WBGene00014450 0 0 0 0
## WBGene00014451 0 0 0 0
## WBGene00010957 0 0 0 0
## WBGene00010958 0 0 0 0
## WBGene00014452 0 0 0 0
## WBGene00014453 0 0 0 0
## L1_GFPminus_rep3 L3_whole_rep1 L3_cells_rep1 L3_GFPplus_rep1
## WBGene00014450 0 0 0 0
## WBGene00014451 0 0 0 0
## WBGene00010957 0 0 0 0
## WBGene00010958 0 0 0 0
## WBGene00014452 0 0 0 0
## WBGene00014453 0 0 0 0
## L3_GFPminus_rep1 L3_whole_rep2 L3_cells_rep2 L3_GFPminus_rep2
## WBGene00014450 0 0 0 0
## WBGene00014451 0 0 0 0
## WBGene00010957 0 0 0 0
## WBGene00010958 0 0 0 0
## WBGene00014452 0 0 0 0
## WBGene00014453 0 0 0 0
## L3_GFPplus_rep2 L3_whole_rep3 L3_cells_rep3 L3_GFPplus_rep3
## WBGene00014450 0 0 0 0
## WBGene00014451 0 0 0 0
## WBGene00010957 0 0 0 0
## WBGene00010958 0 0 0 0
## WBGene00014452 0 0 0 0
## WBGene00014453 0 0 0 0
## L3_GFPminus_rep3
## WBGene00014450 0
## WBGene00014451 0
## WBGene00010957 0
## WBGene00010958 0
## WBGene00014452 0
## WBGene00014453 0# Reorganize the metadata table so the names2 column are now headers
rownames(metadata1)<- metadata1$names
coldata <- metadata1[,c("names", "stage", "sample", "rep")]
rownames(coldata) <- as.vector(metadata1$names)
# make grouping variable
coldata$group <- factor(paste0(coldata$stage, coldata$sample))
coldata## names stage sample rep group
## embryo_cells_rep1 embryo_cells_rep1 embryo cells rep1 embryocells
## embryo_GFPplus_rep1 embryo_GFPplus_rep1 embryo GFPplus rep1 embryoGFPplus
## embryo_GFPminus_rep1 embryo_GFPminus_rep1 embryo GFPminus rep1 embryoGFPminus
## embryo_whole_rep2 embryo_whole_rep2 embryo whole rep2 embryowhole
## embryo_cells_rep2 embryo_cells_rep2 embryo cells rep2 embryocells
## embryo_GFPplus_rep2 embryo_GFPplus_rep2 embryo GFPplus rep2 embryoGFPplus
## embryo_GFPminus_rep2 embryo_GFPminus_rep2 embryo GFPminus rep2 embryoGFPminus
## embryo_whole_rep3 embryo_whole_rep3 embryo whole rep3 embryowhole
## embryo_GFPplus_rep3 embryo_GFPplus_rep3 embryo GFPplus rep3 embryoGFPplus
## embryo_GFPminus_rep3 embryo_GFPminus_rep3 embryo GFPminus rep3 embryoGFPminus
## L1_whole_rep1 L1_whole_rep1 L1 whole rep1 L1whole
## L1_cells_rep1 L1_cells_rep1 L1 cells rep1 L1cells
## L1_GFPplus_rep1 L1_GFPplus_rep1 L1 GFPplus rep1 L1GFPplus
## L1_GFPminus_rep1 L1_GFPminus_rep1 L1 GFPminus rep1 L1GFPminus
## L1_whole_rep2 L1_whole_rep2 L1 whole rep2 L1whole
## L1_cells_rep2 L1_cells_rep2 L1 cells rep2 L1cells
## L1_GFPplus_rep2 L1_GFPplus_rep2 L1 GFPplus rep2 L1GFPplus
## L1_GFPminus_rep2 L1_GFPminus_rep2 L1 GFPminus rep2 L1GFPminus
## L1_whole_rep3 L1_whole_rep3 L1 whole rep3 L1whole
## L1_cells_rep3 L1_cells_rep3 L1 cells rep3 L1cells
## L1_GFPplus_rep3 L1_GFPplus_rep3 L1 GFPplus rep3 L1GFPplus
## L1_GFPminus_rep3 L1_GFPminus_rep3 L1 GFPminus rep3 L1GFPminus
## L3_whole_rep1 L3_whole_rep1 L3 whole rep1 L3whole
## L3_cells_rep1 L3_cells_rep1 L3 cells rep1 L3cells
## L3_GFPplus_rep1 L3_GFPplus_rep1 L3 GFPplus rep1 L3GFPplus
## L3_GFPminus_rep1 L3_GFPminus_rep1 L3 GFPminus rep1 L3GFPminus
## L3_whole_rep2 L3_whole_rep2 L3 whole rep2 L3whole
## L3_cells_rep2 L3_cells_rep2 L3 cells rep2 L3cells
## L3_GFPminus_rep2 L3_GFPminus_rep2 L3 GFPminus rep2 L3GFPminus
## L3_GFPplus_rep2 L3_GFPplus_rep2 L3 GFPplus rep2 L3GFPplus
## L3_whole_rep3 L3_whole_rep3 L3 whole rep3 L3whole
## L3_cells_rep3 L3_cells_rep3 L3 cells rep3 L3cells
## L3_GFPplus_rep3 L3_GFPplus_rep3 L3 GFPplus rep3 L3GFPplus
## L3_GFPminus_rep3 L3_GFPminus_rep3 L3 GFPminus rep3 L3GFPminus# Check that the names match --> Should be TRUE
all(rownames(coldata) == colnames(cts))## [1] TRUEMake DESeqDataSet
Generate the DESeqDataSet. The variables in this design formula will be the type of sample, and the preparation date. This should reduce the variability between the samples based on when they were made.
From the vignette: “In order to benefit from the default settings of the package, you should put the variable of interest at the end of the formula and make sure the control level is the first level.”
The variable of interest is the sample type.
Using DESeqDataSetFromMatrix since I used the program featureCounts.
dds <- DESeqDataSetFromMatrix(countData = cts,
colData = coldata,
design = ~ group)Filter genes with sum counts per million >= 10 across all samples
Visualize read count distribution
raw_count_threshold <- 10
hist(log10(rowSums(counts(dds))), breaks = 50)
abline(v = log10(raw_count_threshold), col = "red", lty = 2)
cpm <- apply(counts(dds),2, function(x) (x/sum(x))*1000000)
hist(log10(rowSums(cpm)))
abline(v = log10(raw_count_threshold), col = "red", lty = 2)
Filter to remove genes with low read counts
keep <- rowSums(cpm) >= raw_count_threshold
dds <- dds[keep,]
dds## class: DESeqDataSet
## dim: 16762 34
## metadata(1): version
## assays(1): counts
## rownames(16762): WBGene00021406 WBGene00021407 ... WBGene00199694
## WBGene00044951
## rowData names(0):
## colnames(34): embryo_cells_rep1 embryo_GFPplus_rep1 ... L3_GFPplus_rep3
## L3_GFPminus_rep3
## colData names(5): names stage sample rep groupFilter to select for protein-coding genes
transcript_type <- read_csv(file = "../01_input/biomaRt_elegans_transcript_biotype.csv")## Rows: 59897 Columns: 4
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (4): Gene stable ID, Genome project, Gene name, Transcript biotype
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.colnames(transcript_type) <- c("WBGeneID", "genome_id", "gene_name", "biotype")
dds <- dds[rownames(dds) %in%
(transcript_type %>%
filter(biotype == "protein_coding") %>%
pull(WBGeneID)),]Perform Differential Expression
dds <- DESeq(dds)## estimating size factors## estimating dispersions## gene-wise dispersion estimates## mean-dispersion relationship## final dispersion estimates## fitting model and testingresultsNames(dds)## [1] "Intercept" "group_embryoGFPminus_vs_embryocells"
## [3] "group_embryoGFPplus_vs_embryocells" "group_embryowhole_vs_embryocells"
## [5] "group_L1cells_vs_embryocells" "group_L1GFPminus_vs_embryocells"
## [7] "group_L1GFPplus_vs_embryocells" "group_L1whole_vs_embryocells"
## [9] "group_L3cells_vs_embryocells" "group_L3GFPminus_vs_embryocells"
## [11] "group_L3GFPplus_vs_embryocells" "group_L3whole_vs_embryocells"Sample-to-sample distance matrix
vsd <- vst(dds, blind = FALSE)
sampleDists <- dist(t(assay(vsd)))
sampleDistMatrix <- as.matrix(sampleDists)
rownames(sampleDistMatrix) <- vsd$names
colnames(sampleDistMatrix) <- NULL
colors <- colorRampPalette( rev(brewer.pal(9, "Blues")) )(255)
mega_cor_plot <- pheatmap(sampleDistMatrix,
clustering_distance_rows = sampleDists,
clustering_distance_cols = sampleDists,
col = colors)
mega_cor_plot
myPDFplot(plot = mega_cor_plot, name = "FACS_Correlation_All_Samples", height = 4.5, width = 6, plotdir = "../03_output/plots/Correlation_Matrix/")## quartz_off_screen
## 2Figure S2
all_sorted_samples_cor <- vsd.corr.per.stage("GFPplus|GFPminus", "Correlation of FACS isolated GFP+ and GFP- samples")
myPDFplot(plot = all_sorted_samples_cor, name = "All_Stage_FACS_Correlation_Sorted_Samples", height = 4, width = 6, plotdir = "../03_output/plots/Correlation_Matrix/")## quartz_off_screen
## 2Remove L1 rep 2
remove_samples <- c("L1_whole_rep2", "L1_cells_rep2", "L1_GFPplus_rep2", "L1_GFPminus_rep2")
coldata <- coldata %>% filter(!names %in% remove_samples)
dds <- dds[,!colnames(dds)%in% remove_samples]Perform Differential Expression
dds <- DESeq(dds)## using pre-existing size factors## estimating dispersions## found already estimated dispersions, replacing these## gene-wise dispersion estimates## mean-dispersion relationship## final dispersion estimates## fitting model and testingresultsNames(dds)## [1] "Intercept" "group_embryoGFPminus_vs_embryocells"
## [3] "group_embryoGFPplus_vs_embryocells" "group_embryowhole_vs_embryocells"
## [5] "group_L1cells_vs_embryocells" "group_L1GFPminus_vs_embryocells"
## [7] "group_L1GFPplus_vs_embryocells" "group_L1whole_vs_embryocells"
## [9] "group_L3cells_vs_embryocells" "group_L3GFPminus_vs_embryocells"
## [11] "group_L3GFPplus_vs_embryocells" "group_L3whole_vs_embryocells"Counts table output for supplemental datatable
# raw counts
write.table(as.data.frame(counts(dds, normalized = FALSE)) %>% rownames_to_column(var = "WBGeneID"),
file = "../03_output/count_tables_for_sup/intestine_FACS_RNAseq_raw_counts.txt",
sep = "\t", quote = FALSE, row.names = FALSE)
# normalized counts
write.table(as.data.frame(counts(dds, normalized = TRUE)) %>% rownames_to_column(var = "WBGeneID"),
file = "../03_output/count_tables_for_sup/intestine_FACS_RNAseq_norm_counts.txt",
sep = "\t", quote = FALSE, row.names = FALSE)
# rlog transformed counts
write.table(as.data.frame(assay(rlog(dds, blind=FALSE))) %>% rownames_to_column(var = "WBGeneID"),
file = "../03_output/count_tables_for_sup/intestine_FACS_RNAseq_rlog_counts.txt",
sep = "\t", quote = FALSE, row.names = FALSE)## rlog() may take a few minutes with 30 or more samples,
## vst() is a much faster transformationWithin-stage pairwise comparisons
Set cutoff values
thresh = 1
sig = 0.01Embryo stage pairwise comparisons
embryo_alt_hyp_res_df <- alt_hyp_res_df("embryo", thresh = thresh, sig = sig)de_category_MA_plot(embryo_alt_hyp_res_df, paste("Embryo differentially expressed genes\nlfc = ",thresh," & padj < ",sig, sep = ""))## Warning: Removed 91 rows containing missing values (geom_point).## Warning: Removed 32 rows containing missing values (geom_point).
de_category_bar_plot(embryo_alt_hyp_res_df, paste("Embryo differentially expressed genes\nlfc = ",thresh," & padj < ",sig, sep = ""))## `summarise()` has grouped output by 'label'. You can override using the
## `.groups` argument.
## L1 stage pairwise comparisons
L1_alt_hyp_res_df<- alt_hyp_res_df("L1", thresh = thresh, sig = sig)de_category_MA_plot(L1_alt_hyp_res_df, paste("L1 differentially expressed genes\nlfc = ",thresh," & padj < ",sig, sep = ""))## Warning: Removed 64 rows containing missing values (geom_point).## Warning: Removed 31 rows containing missing values (geom_point).
de_category_bar_plot(L1_alt_hyp_res_df, paste("L1 differentially expressed genes\nlfc = ",thresh," & padj < ",sig, sep = ""))## `summarise()` has grouped output by 'label'. You can override using the
## `.groups` argument.
## L3 stage pairwise comparisons
L3_alt_hyp_res_df<- alt_hyp_res_df("L3", thresh = thresh, sig = sig)de_category_MA_plot(L3_alt_hyp_res_df, paste("L3 differentially expressed genes\nlfc = ",thresh," & padj < ",sig, sep = ""))
de_category_bar_plot(L3_alt_hyp_res_df, paste("L3 differentially expressed genes\nlfc = ",thresh," & padj < ",sig, sep = ""))## `summarise()` has grouped output by 'label'. You can override using the
## `.groups` argument.
# Figure S3
Whole worm vs. dissociated cells analysis
cells_vs_whole_bar_df <- L1_alt_hyp_res_df %>% filter(isDE == TRUE, label == "whole_vs_cells") %>% group_by(type) %>% summarize(genes = n()) %>% mutate(stage = "L1") %>%
bind_rows(
L3_alt_hyp_res_df %>% filter(isDE == TRUE, label == "whole_vs_cells") %>% group_by(type) %>% summarize(genes = n()) %>% mutate(stage = "L3")
) %>%
bind_rows(
embryo_alt_hyp_res_df %>% filter(isDE == TRUE, label == "whole_vs_cells") %>% group_by(type) %>% summarize(genes = n()) %>% mutate(stage = "embryo")
)
cells_vs_whole_bar_df## # A tibble: 9 × 3
## type genes stage
## <chr> <int> <chr>
## 1 greater 494 L1
## 2 less 11 L1
## 3 lessAbs 1117 L1
## 4 greater 419 L3
## 5 less 9 L3
## 6 lessAbs 2058 L3
## 7 greater 7 embryo
## 8 less 1 embryo
## 9 lessAbs 1187 embryocells_vs_whole_bar_plot <- cells_vs_whole_bar_df %>%
ggplot(aes(x = type, y = genes, fill = stage, label = genes)) +
geom_bar(stat = "identity", position = position_dodge(width = 0.7), width = 0.5, color = "black") +
geom_text(hjust = -0.25, position = position_dodge(width = 0.7)) +
scale_fill_manual(values = c("black", "grey", "white")) +
theme_classic() +
coord_flip()
cells_vs_whole_bar_plot
ggsave(plot = cells_vs_whole_bar_plot, filename = "../03_output/plots/cells_vs_whole/cells_vs_whole_bar_plot.pdf", width = 5.5, height = 3.5) cells_vs_whole_MA_plot <- de_category_MA_plot(
L1_alt_hyp_res_df %>% filter(label == "whole_vs_cells") %>% mutate(label = comparison) %>%
bind_rows(
L3_alt_hyp_res_df %>% filter(label == "whole_vs_cells") %>% mutate(label = comparison)
) %>%
bind_rows(
embryo_alt_hyp_res_df %>% filter(label == "whole_vs_cells") %>% mutate(label = comparison)
),
title = NULL
)
cells_vs_whole_MA_plot## Warning: Removed 10 rows containing missing values (geom_point).## Warning: Removed 2 rows containing missing values (geom_point).
ggsave(plot = cells_vs_whole_MA_plot, filename = "../03_output/plots/cells_vs_whole/cells_vs_whole_MA_plot.pdf", width = 6, height = 3)## Warning: Removed 10 rows containing missing values (geom_point).
## Removed 2 rows containing missing values (geom_point).DESeq2 builtin plotting function
res_embryoGFPplus_vs_embryoGFPminus <- results(dds, contrast = c("group", "embryoGFPplus", "embryoGFPminus"))
res_L1GFPplus_vs_L1_GFPminus <- results(dds, contrast = c("group", "L1GFPplus", "L1GFPminus"))
res_L3GFPplus_vs_L3_GFPminus <- results(dds, contrast = c("group", "L3GFPplus", "L3GFPminus"))res_embryoGFPplus_vs_embryoGFPminus_ashr <- lfcShrink(dds, contrast = c("group", "embryoGFPplus", "embryoGFPminus"), type = "ashr")## using 'ashr' for LFC shrinkage. If used in published research, please cite:
## Stephens, M. (2016) False discovery rates: a new deal. Biostatistics, 18:2.
## https://doi.org/10.1093/biostatistics/kxw041res_L1GFPplus_vs_L1GFPminus_ashr <- lfcShrink(dds, contrast = c("group", "L1GFPplus", "L1GFPminus"), type = "ashr")## using 'ashr' for LFC shrinkage. If used in published research, please cite:
## Stephens, M. (2016) False discovery rates: a new deal. Biostatistics, 18:2.
## https://doi.org/10.1093/biostatistics/kxw041res_L3GFPplus_vs_L3GFPminus_ashr <- lfcShrink(dds, contrast = c("group", "L3GFPplus", "L3GFPminus"), type = "ashr")## using 'ashr' for LFC shrinkage. If used in published research, please cite:
## Stephens, M. (2016) False discovery rates: a new deal. Biostatistics, 18:2.
## https://doi.org/10.1093/biostatistics/kxw041Export rlog counts
all_samples_rld <- rlog(dds)## rlog() may take a few minutes with 30 or more samples,
## vst() is a much faster transformationwrite_rds(all_samples_rld, file = "../03_output/rlog_counts/all_samples_rlog_counts.rds")
write_tsv(as.data.frame(assay(all_samples_rld)) %>% rownames_to_column(var = "WBGeneID"), file = "../03_output/rlog_counts/all_samples_rlog_counts.tsv")
write_tsv(as.data.frame(assay(all_samples_rld)) %>% rownames_to_column(var = "WBGeneID") %>% select(WBGeneID, contains("GFPplus")), file = "../03_output/rlog_counts/GFPplus_samples_rlog_counts.tsv")
all_samples_rld <- read_rds(file = "../03_output/rlog_counts/all_samples_rlog_counts.rds")
all_samples_rld_df <- read_tsv(file = "../03_output/rlog_counts/all_samples_rlog_counts.tsv")## Rows: 15627 Columns: 31## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## chr (1): WBGeneID
## dbl (30): embryo_cells_rep1, embryo_GFPplus_rep1, embryo_GFPminus_rep1, embr...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.head(all_samples_rld_df)## # A tibble: 6 × 31
## WBGeneID embryo_cells_re… embryo_GFPplus_… embryo_GFPminus… embryo_whole_re…
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 WBGene000… 6.21 6.19 6.70 5.52
## 2 WBGene000… 1.94 3.43 2.92 2.99
## 3 WBGene000… 3.63 3.35 0.907 1.30
## 4 WBGene000… 1.15 3.27 0.760 1.95
## 5 WBGene000… 9.82 10.1 9.84 8.95
## 6 WBGene000… 1.59 1.06 1.48 2.14
## # … with 26 more variables: embryo_cells_rep2 <dbl>, embryo_GFPplus_rep2 <dbl>,
## # embryo_GFPminus_rep2 <dbl>, embryo_whole_rep3 <dbl>,
## # embryo_GFPplus_rep3 <dbl>, embryo_GFPminus_rep3 <dbl>, L1_whole_rep1 <dbl>,
## # L1_cells_rep1 <dbl>, L1_GFPplus_rep1 <dbl>, L1_GFPminus_rep1 <dbl>,
## # L1_whole_rep3 <dbl>, L1_cells_rep3 <dbl>, L1_GFPplus_rep3 <dbl>,
## # L1_GFPminus_rep3 <dbl>, L3_whole_rep1 <dbl>, L3_cells_rep1 <dbl>,
## # L3_GFPplus_rep1 <dbl>, L3_GFPminus_rep1 <dbl>, L3_whole_rep2 <dbl>, …Export intestine enrichment assignment
embryo_intestine_gene_cats <- embryo_alt_hyp_res_df %>%
drop_na(padj) %>%
filter(label == "GFPplus_vs_GFPminus", padj < sig) %>%
select(WBGeneID, altHyp = "type") %>%
mutate(intestine_expression = case_when(
altHyp == "greater" ~ "enriched",
altHyp == "less" ~ "depleted",
altHyp == "lessAbs" ~ "equal")) %>%
mutate(intestine_expression = fct_relevel(intestine_expression, c("enriched", "equal", "depleted")))
L1_intestine_gene_cats <- L1_alt_hyp_res_df %>%
filter(label == "GFPplus_vs_GFPminus", padj < sig) %>%
select(WBGeneID, altHyp = "type") %>%
mutate(intestine_expression = case_when(
altHyp == "greater" ~ "enriched",
altHyp == "less" ~ "depleted",
altHyp == "lessAbs" ~ "equal"))%>%
mutate(intestine_expression = fct_relevel(intestine_expression, c("enriched", "equal", "depleted")))
L3_intestine_gene_cats <- L3_alt_hyp_res_df %>%
filter(label == "GFPplus_vs_GFPminus", padj < sig) %>%
select(WBGeneID, altHyp = "type") %>%
mutate(intestine_expression = case_when(
altHyp == "greater" ~ "enriched",
altHyp == "less" ~ "depleted",
altHyp == "lessAbs" ~ "equal"))%>%
mutate(intestine_expression = fct_relevel(intestine_expression, c("enriched", "equal", "depleted")))
write_csv(x = embryo_intestine_gene_cats, file = "../03_output/intestine_gene_categories/embryo_intestine_gene_categories.csv")
write_csv(x = L1_intestine_gene_cats, file = "../03_output/intestine_gene_categories/L1_intestine_gene_categories.csv")
write_csv(x = L3_intestine_gene_cats, file = "../03_output/intestine_gene_categories/L3_intestine_gene_categories.csv")Average embryo GFP+ sample reads
Make function
thresh = 1
sig = 0.01
rlog_status <- function(stage, res, hyp_df){
all_samples_rld_df %>%
select(WBGeneID, contains(paste(stage,"GFPplus", sep = "_"))) %>%
pivot_longer(cols = contains("GFPplus"), values_to = "rlog_counts") %>%
separate(name, sep = "_", into = c("stage", "sample", "rep")) %>%
group_by(WBGeneID) %>%
summarise(mean.rlog.counts = mean(rlog_counts), var.rlog.counts = var(rlog_counts)) %>%
left_join(hyp_df %>% filter(label == "GFPplus_vs_GFPminus") %>% select(WBGeneID, type, isDE), by = "WBGeneID")
}embryo_rlog_status_df <- rlog_status(stage = "embryo", res = res_embryoGFPplus_vs_embryoGFPminus, hyp_df = embryo_alt_hyp_res_df)
head(embryo_rlog_status_df)## # A tibble: 6 × 5
## WBGeneID mean.rlog.counts var.rlog.counts type isDE
## <chr> <dbl> <dbl> <chr> <lgl>
## 1 WBGene00000001 10.1 0.00518 greater FALSE
## 2 WBGene00000001 10.1 0.00518 less FALSE
## 3 WBGene00000001 10.1 0.00518 lessAbs FALSE
## 4 WBGene00000002 10.9 0.135 greater FALSE
## 5 WBGene00000002 10.9 0.135 less FALSE
## 6 WBGene00000002 10.9 0.135 lessAbs FALSEwrite_csv(embryo_rlog_status_df, file = "../03_output/embryo_GFPplus_rlog_counts_status_df.csv", col_names = TRUE)
embryo_rlog_status_df <- read_csv(file = "../03_output/embryo_GFPplus_rlog_counts_status_df.csv", col_names = TRUE)## Rows: 40203 Columns: 5
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (2): WBGeneID, type
## dbl (2): mean.rlog.counts, var.rlog.counts
## lgl (1): isDE
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.L1_rlog_status_df <- rlog_status(stage = "L1", res = res_L1GFPplus_vs_L1_GFPminus, hyp_df = L1_alt_hyp_res_df)
write_csv(L1_rlog_status_df, file = "../03_output/L1_GFPplus_rlog_counts_status_df.csv", col_names = TRUE)
L1_rlog_status_df <- read_csv(file = "../03_output/L1_GFPplus_rlog_counts_status_df.csv", col_names = TRUE)## Rows: 38690 Columns: 5
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (2): WBGeneID, type
## dbl (2): mean.rlog.counts, var.rlog.counts
## lgl (1): isDE
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.head(L1_rlog_status_df)## # A tibble: 6 × 5
## WBGeneID mean.rlog.counts var.rlog.counts type isDE
## <chr> <dbl> <dbl> <chr> <lgl>
## 1 WBGene00000001 10.1 0.0190 greater FALSE
## 2 WBGene00000001 10.1 0.0190 less FALSE
## 3 WBGene00000001 10.1 0.0190 lessAbs TRUE
## 4 WBGene00000002 8.34 0.142 greater FALSE
## 5 WBGene00000002 8.34 0.142 less FALSE
## 6 WBGene00000002 8.34 0.142 lessAbs FALSEL3_rlog_status_df <- rlog_status(stage = "L3", res = res_L3GFPplus_vs_L3_GFPminus, hyp_df = L3_alt_hyp_res_df)
write_csv(L3_rlog_status_df, file = "../03_output/L3_GFPplus_rlog_counts_status_df.csv", col_names = TRUE)
L3_rlog_status_df <- read_csv(file = "../03_output/L3_GFPplus_rlog_counts_status_df.csv", col_names = TRUE)## Rows: 39901 Columns: 5
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (2): WBGeneID, type
## dbl (2): mean.rlog.counts, var.rlog.counts
## lgl (1): isDE
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.head(L3_rlog_status_df)## # A tibble: 6 × 5
## WBGeneID mean.rlog.counts var.rlog.counts type isDE
## <chr> <dbl> <dbl> <chr> <lgl>
## 1 WBGene00000001 10.1 0.0919 greater FALSE
## 2 WBGene00000001 10.1 0.0919 less FALSE
## 3 WBGene00000001 10.1 0.0919 lessAbs TRUE
## 4 WBGene00000002 8.99 0.352 greater FALSE
## 5 WBGene00000002 8.99 0.352 less FALSE
## 6 WBGene00000002 8.99 0.352 lessAbs FALSEIntesinte expression across development
UpSet Plot
intestine_enriched_genes <- data.frame(embryo_rlog_status_df, stage = "embryo") %>% bind_rows(data.frame(L1_rlog_status_df, stage = "L1"), data.frame(L3_rlog_status_df, stage = "L3")) %>%filter(type == "greater", isDE == TRUE) %>% select(WBGeneID, stage)
stage_list<- list(embryo = filter(intestine_enriched_genes, stage == "embryo")$WBGeneID,
L1 = filter(intestine_enriched_genes, stage == "L1")$WBGeneID,
L3 = filter(intestine_enriched_genes, stage == "L3")$WBGeneID)
comb_mat <-make_comb_mat(stage_list)
UpSet(comb_mat)
comb_size(comb_mat)## 111 110 101 011 100 010 001
## 848 114 130 302 644 213 154Data output
# per stage GFP+ vs GFP- differential expression
write_csv(res_to_df(res_embryoGFPplus_vs_embryoGFPminus), file = "../03_output/pairwise_DE_results/res_embryoGFPplus_vs_embryoGFPminus.csv")
write_csv(res_to_df(res_L1GFPplus_vs_L1_GFPminus), file = "../03_output/pairwise_DE_results/res_L1GFPplus_vs_L1GFPminus.csv")
write_csv(res_to_df(res_L3GFPplus_vs_L3_GFPminus), file = "../03_output/pairwise_DE_results/res_L3GFPplus_vs_L3GFPminus.csv")
# per stage GFP+ vs GFP- differential expression with log2FC shrink (visualization, ranking)
write_csv(res_to_df(res_embryoGFPplus_vs_embryoGFPminus_ashr), file = "../03_output/pairwise_shrunk_DE_results/res_embryoGFPplus_vs_embryoGFPminus_ashr_shrunk.csv")
write_csv(res_to_df(res_L1GFPplus_vs_L1GFPminus_ashr), file = "../03_output/pairwise_shrunk_DE_results/res_L1GFPplus_vs_L1GFPminus_ashr_shrunk.csv")
write_csv(res_to_df(res_L3GFPplus_vs_L3GFPminus_ashr), file = "../03_output/pairwise_shrunk_DE_results/res_L3GFPplus_vs_L3GFPminus_ashr_shrunk.csv")
# per stage enrichment annotation
write_csv(embryo_alt_hyp_res_df %>% filter(label == "GFPplus_vs_GFPminus") %>% select(WBGeneID:type, isDE), file = "../03_output/embryo_GFPplus_vs_GFPminus_alt_hyp_res.csv")
write_csv(L1_alt_hyp_res_df %>% filter(label == "GFPplus_vs_GFPminus") %>% select(WBGeneID:type, isDE), file = "../03_output/L1_GFPplus_vs_GFPminus_alt_hyp_res.csv")
write_csv(L3_alt_hyp_res_df %>% filter(label == "GFPplus_vs_GFPminus") %>% select(WBGeneID:type, isDE), file = "../03_output/L3_GFPplus_vs_GFPminus_alt_hyp_res.csv")Plot output
res_embryoGFPplus_vs_embryoGFPminus_ashr <- read_csv(file = "../03_output/pairwise_shrunk_DE_results/res_embryoGFPplus_vs_embryoGFPminus_ashr_shrunk.csv")## Rows: 15627 Columns: 6
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): WBGeneID
## dbl (5): baseMean, log2FoldChange, lfcSE, pvalue, padj
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.embryo_intestine_gene_cats <- read_csv(file = "../03_output/intestine_gene_categories/embryo_intestine_gene_categories.csv")## Rows: 3142 Columns: 3
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (3): WBGeneID, altHyp, intestine_expression
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.res_embryoGFPplus_vs_embryoGFPminus_ashr %>% left_join(embryo_intestine_gene_cats, by = "WBGeneID") %>% drop_na(intestine_expression) %>% filter(intestine_expression == "enriched") %>%
ggplot(aes(x = log10(baseMean), y =log2FoldChange)) +
geom_point(data = res_embryoGFPplus_vs_embryoGFPminus_ashr %>% filter(!(WBGeneID %in% embryo_intestine_gene_cats$WBGeneID)),
shape = 16, alpha = 0.5, stroke = 0, size = 1, color = "grey") +
geom_point(shape = 16, alpha = 0.5, stroke = 0, size = 1, aes(color = intestine_expression)) +
theme_bw() +
labs(x = "log10(mean RNA abundance)",
y = "log2(embryo GFP+ reads/embryo GFP- reads)")
ggsave(filename = "../03_output/plots/Intestine_Expression_Category/L3_GFPplus_vs_GFPminus_shrunk_MAplot_enriched_only.jpg", width = 5, height = 3)
res_to_df(res_embryoGFPplus_vs_embryoGFPminus_ashr) %>% left_join(embryo_intestine_gene_cats, by = "WBGeneID") %>% drop_na(intestine_expression) %>% group_by(intestine_expression) %>% summarise(genes = n()) %>% ggplot(aes(x = intestine_expression, y = genes, fill = intestine_expression, label = genes)) +
geom_bar(stat = "identity") +
geom_text(vjust = -0.25) +
theme_bw()
single_MA_plot <- function(in_res, in_cats){
res_to_df(in_res) %>% left_join(in_cats, by = "WBGeneID") %>% drop_na(intestine_expression) %>%
filter(log2FoldChange > -20, intestine_expression == "enriched") %>%
ggplot(aes(x = log10(baseMean), y =log2FoldChange)) +
geom_point(data = res_to_df(in_res) %>% filter(!(WBGeneID %in% in_cats$WBGeneID)),
shape = 16, alpha = 0.5, stroke = 0, size = 1, color = "grey") +
geom_point(shape = 16, alpha = 0.5, stroke = 0, size = 1, aes(color = intestine_expression)) +
theme_bw() +
xlim(0.5,5.5) +
ylim(-6, 10) +
labs(x = "log10(mean RNA abundance)",
y = "log2(GFP+ reads/GFP- reads)",
title = paste("data: ", deparse(substitute(in_res)), sep = ""))
}
single_cat_bar <- function(in_res, in_cats){
res_to_df(in_res) %>% left_join(in_cats, by = "WBGeneID") %>% drop_na(intestine_expression) %>% group_by(intestine_expression) %>% summarise(genes = n()) %>% ggplot(aes(x = intestine_expression, y = genes, fill = intestine_expression, label = genes)) +
geom_bar(stat = "identity") +
geom_text(vjust = -0.25) +
ggtitle(paste("data: ", deparse(substitute(in_cats)), sep = "")) +
theme_bw()
}embryo_shrunk_MA <- single_MA_plot(res_embryoGFPplus_vs_embryoGFPminus_ashr, embryo_intestine_gene_cats)
embryo_shrunk_MA## Warning: Removed 3 rows containing missing values (geom_point).## Warning: Removed 10 rows containing missing values (geom_point).
ggsave(plot = embryo_shrunk_MA, file = "../03_output/plots/Intestine_Expression_Category/embryo_GFPplus_vs_GFPminus_shrunk_MAplot.pdf", width = 5, height = 3)## Warning: Removed 3 rows containing missing values (geom_point).
## Removed 10 rows containing missing values (geom_point).embryo_category_bar <- single_cat_bar(res_embryoGFPplus_vs_embryoGFPminus_ashr, embryo_intestine_gene_cats)
embryo_category_bar
ggsave(plot = embryo_category_bar, file = "../03_output/plots/Intestine_Expression_Category/embryo_GFPplus_vs_GFPminus_Expression_Category_barplot.pdf", width = 5, height = 3)L1_shrunk_MA <- single_MA_plot(res_L1GFPplus_vs_L1GFPminus_ashr, L1_intestine_gene_cats)
L1_shrunk_MA## Warning: Removed 3 rows containing missing values (geom_point).
ggsave(plot = L1_shrunk_MA, file = "../03_output/plots/Intestine_Expression_Category/L1_GFPplus_vs_GFPminus_shrunk_MAplot.pdf", width = 5, height = 3)## Warning: Removed 3 rows containing missing values (geom_point).L1_category_bar <- single_cat_bar(res_L1GFPplus_vs_L1GFPminus_ashr, L1_intestine_gene_cats)
L1_category_bar
ggsave(plot = L1_category_bar, file = "../03_output/plots/Intestine_Expression_Category/L1_GFPplus_vs_GFPminus_Expression_Category_barplot.pdf", width = 5, height = 3)L3_shrunk_MA <- single_MA_plot(res_L3GFPplus_vs_L3GFPminus_ashr, L3_intestine_gene_cats)
L3_shrunk_MA## Warning: Removed 3 rows containing missing values (geom_point).
ggsave(plot = L3_shrunk_MA, file = "../03_output/plots/Intestine_Expression_Category/L3_GFPplus_vs_GFPminus_shrunk_MAplot.pdf", width = 5, height = 3)## Warning: Removed 3 rows containing missing values (geom_point).L3_category_bar <- single_cat_bar(res_L3GFPplus_vs_L3GFPminus_ashr, L3_intestine_gene_cats)
L3_category_bar
ggsave(plot = L3_category_bar, file = "../03_output/plots/Intestine_Expression_Category/L3_GFPplus_vs_GFPminus_Expression_Category_barplot.pdf", width = 5, height = 3)Alternative Hypothesis output
thresh <- 1
sig <- 0.01
res_embryoGFP_alHyp_greater <- results(dds, contrast = c("group", "embryoGFPplus", "embryoGFPminus"), lfcThreshold=thresh, altHypothesis = "greater", alpha = sig)
write_csv(x = res_to_df(res_embryoGFP_alHyp_greater), file = "../03_output/res_embryoGFP_alHyp_greater.csv")
res_L1GFP_alHyp_greater <- results(dds, contrast = c("group", "L1GFPplus", "L1GFPminus"), lfcThreshold=thresh, altHypothesis = "greater", alpha = sig)
write_csv(x = res_to_df(res_L1GFP_alHyp_greater), file = "../03_output/res_L1GFP_alHyp_greater.csv")
res_L3GFP_alHyp_greater <- results(dds, contrast = c("group", "L3GFPplus", "L3GFPminus"), lfcThreshold=thresh, altHypothesis = "greater", alpha = sig)
write_csv(x = res_to_df(res_L3GFP_alHyp_greater), file = "../03_output/res_L3GFP_alHyp_greater.csv")Session info
sessionInfo()## R version 4.1.0 (2021-05-18)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Catalina 10.15.7
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] grid parallel stats4 stats graphics grDevices utils
## [8] datasets methods base
##
## other attached packages:
## [1] ComplexHeatmap_2.8.0 InterMineR_1.14.1
## [3] ashr_2.2-54 apeglm_1.14.0
## [5] forcats_0.5.1 stringr_1.4.0
## [7] dplyr_1.0.8 purrr_0.3.4
## [9] readr_2.1.2 tidyr_1.2.0
## [11] tibble_3.1.6 ggplot2_3.3.5
## [13] tidyverse_1.3.1 pheatmap_1.0.12
## [15] RColorBrewer_1.1-3 corrplot_0.92
## [17] DESeq2_1.32.0 SummarizedExperiment_1.22.0
## [19] Biobase_2.52.0 MatrixGenerics_1.4.3
## [21] matrixStats_0.61.0 GenomicRanges_1.44.0
## [23] GenomeInfoDb_1.28.4 IRanges_2.26.0
## [25] S4Vectors_0.30.2 BiocGenerics_0.38.0
##
## loaded via a namespace (and not attached):
## [1] readxl_1.4.0 backports_1.4.1 circlize_0.4.14
## [4] plyr_1.8.7 igraph_1.3.0 splines_4.1.0
## [7] BiocParallel_1.26.2 digest_0.6.29 invgamma_1.1
## [10] foreach_1.5.2 htmltools_0.5.2 SQUAREM_2021.1
## [13] fansi_1.0.3 magrittr_2.0.3 memoise_2.0.1
## [16] cluster_2.1.3 doParallel_1.0.17 tzdb_0.3.0
## [19] Biostrings_2.60.2 annotate_1.70.0 modelr_0.1.8
## [22] vroom_1.5.7 bdsmatrix_1.3-4 colorspace_2.0-3
## [25] blob_1.2.3 rvest_1.0.2 haven_2.4.3
## [28] xfun_0.30 crayon_1.5.1 RCurl_1.98-1.6
## [31] jsonlite_1.8.0 genefilter_1.74.1 survival_3.3-1
## [34] iterators_1.0.14 glue_1.6.2 gtable_0.3.0
## [37] zlibbioc_1.38.0 XVector_0.32.0 GetoptLong_1.0.5
## [40] DelayedArray_0.18.0 shape_1.4.6 scales_1.2.0
## [43] mvtnorm_1.1-3 DBI_1.1.2 Rcpp_1.0.8.3
## [46] xtable_1.8-4 emdbook_1.3.12 clue_0.3-60
## [49] bit_4.0.4 sqldf_0.4-11 truncnorm_1.0-8
## [52] httr_1.4.2 ellipsis_0.3.2 farver_2.1.0
## [55] pkgconfig_2.0.3 XML_3.99-0.9 sass_0.4.1
## [58] dbplyr_2.1.1 locfit_1.5-9.5 utf8_1.2.2
## [61] RJSONIO_1.3-1.6 labeling_0.4.2 tidyselect_1.1.2
## [64] rlang_1.0.2 AnnotationDbi_1.54.1 munsell_0.5.0
## [67] cellranger_1.1.0 tools_4.1.0 cachem_1.0.6
## [70] cli_3.2.0 gsubfn_0.7 generics_0.1.2
## [73] RSQLite_2.2.12 broom_0.8.0 evaluate_0.15
## [76] fastmap_1.1.0 yaml_2.3.5 knitr_1.38
## [79] bit64_4.0.5 fs_1.5.2 KEGGREST_1.32.0
## [82] xml2_1.3.3 compiler_4.1.0 rstudioapi_0.13
## [85] png_0.1-7 reprex_2.0.1 geneplotter_1.70.0
## [88] bslib_0.3.1 stringi_1.7.6 highr_0.9
## [91] lattice_0.20-45 Matrix_1.4-1 vctrs_0.4.0
## [94] pillar_1.7.0 lifecycle_1.0.1 jquerylib_0.1.4
## [97] GlobalOptions_0.1.2 bitops_1.0-7 irlba_2.3.5
## [100] R6_2.5.1 codetools_0.2-18 MASS_7.3-56
## [103] assertthat_0.2.1 chron_2.3-56 proto_1.0.0
## [106] rjson_0.2.21 withr_2.5.0 GenomeInfoDbData_1.2.6
## [109] hms_1.1.1 coda_0.19-4 rmarkdown_2.13
## [112] Cairo_1.5-15 mixsqp_0.3-43 bbmle_1.0.24
## [115] numDeriv_2016.8-1.1 lubridate_1.8.0